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Abstract. Viscous heating plays an important role in the dy- 
namics of fluids with strongly temperature-dependent viscos- 
ity because of the coupling between the energy and momen- 
tum equations. The heat generated by viscous friction pro- 
duces a local temperature increase near the tube walls with 
a consequent decrease of the viscosity which may dramat- 
ically change the temperature and velocity profiles. These 
processes are mainly controlled by the Peclet number, the 
Nahme number, the flow rate and the thermal boundary con- 
ditions. The problem of viscous heating in fluids was in- 
vestigated in the past for its practical interest in the polymer 
industry, and was invoked to explain some rheological be- 
haviours of silicate melts, but was not completely applied to 
study magma flows. In this paper we focus on the thermal 
and mechanical effects caused by viscous heating in tubes of 
finite lengths. We find that in magma flows at high Nahme 
number and typical flow rates, viscous heating is responsi- 
ble for the evolution from Poiseuille flow, with a uniform 
temperature distribution at the inlet, to a plug flow with a 
hotter layer near the walls. When the temperature gradients 
induced by viscous heating are very pronounced, local insta- 
bilities may occur and the triggering of secondary flows is 
possible. For completeness, this paper also describes magma 
flow in infinitely long tubes both at steady state and in tran- 
sient phase. 



1 Introduction 

A common feature of lava flows is the presence of a complex 
network of lava tubes, observed both in "pahoehoe" (Peter- 
son et al., 1994) and 'a'a lava fields (Calvari and Pinkerton, 
1998, 1999). The tube system acts as an efficient pathway 
for lava from the main vent to the flow fronts, allowing great 
distances to be reached. This is mostly due to the decrease of 
radiative cooling, as a consequence of a much lower tempera- 
ture of the outer crust of the tubes as compared to the temper- 
ature of the lava surface flowing in open channels (Dragoni, 
1989; Dragoni et al., 1995). 



The effects of heat generation by viscous friction for lava 
flows in tubes and channels were often neglected in previous 
models. In this work we show that these effects can play 
an important role in the dynamics of fluids with strongly 
temperature-dependent viscosity such as silicate melts and 
polymers. In fact, in these fluids, viscous friction generates 
a local increase in temperature near the tube walls with con- 
sequent viscosity decrease and increase of the flow velocity. 
This velocity increase produces a further growth of the lo- 
cal temperature. As we will see later, there are some criti- 
cal values of the parameters that control this process above 
which this feedback cannot converge. In this case the one- 
dimensional laminar solution, valid in the limit of infinitely 
long pipe, cannot exist even for low Reynolds numbers. In 
pipes with finite length, viscous heating governs the evolu- 
tion from Poiseuille regime with a uniform temperature dis- 
tribution at the conduit inlet, to a plug flow with a hotter 
boundary layer near the walls downstream. This effect could 
be observed even in fluids with null yield strength. When the 
temperature gradients, induced by viscous heating are very 
pronounced, local instabilities occur and triggering of sec- 
ondary flows is possible. 

Viscous heating was previously invoked to describe the 
rheological behaviour of basalts (Shaw, 1969) and to explain 
some instabilities in the mantle (Anderson and Perkins, 1974; 
Zhao and Yuen, 1987; Larsen et al., 1995; Hansen and Yuen, 
1996). 

Moreover Fujii and Uyeda (1974), adopting a model only 
valid for infinitely long tubes (by the thermal point of view), 
explained the size of intrusive dikes in volcanoes. 

Because of the low thermal conductivity of silicate melts, 
the temperature field shows a strong radial gradient and vis- 
cosity layering. Flows with layers of different viscosity were 
investigated in the past, also for their practical interest, and 
it is known that they are not always stable (Yih, 1967; Craik, 
1969; Renardy and Joseph, 1985; Renardy, 1987). The insta- 
bilities are generated by the viscosity contrasts, and are sim- 
ilar to the Kelvin-Helmholtz instabilities triggered by den- 
sity gradients. In Couette flows of fluids with temperature- 
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dependent viscosity, Sukanek and Laurence (1974) and Yueh 
and Weng (1996) found one instability mode related to the 
viscosity gradient at low Reynolds numbers. This kind of 
instability shows a local character (Pearson, 1976) and is ex- 
pected to occur also in lava flows. 

As we will see later, these processes are controlled prin- 
cipally by three parameters: the Peclet number (Pe), the 
Nahme number (Na , also called Brinkman number), and 
the non-dimensional flow rate (q): 

Pe = pcUH/k; Na a = p U 2 f3/k; q = p Q/(pgH 3 )(l) 

with p density, c specific heat, U mean velocity, H tube ra- 
dius, k thermal conductivity, p,Q reference viscosity (Nao is 
based on this value), (3 rheological parameter (see eq.5) and 
Q flow rate per unit length (Q = U H). 

Viscous heating in lava and magma flows is responsible 
for effects not described by simple isothermal models, and 
may help in the understanding of some phenomena observed 
in lava flows, which are not yet clearly described. These in- 
clude the surprising temperature increase at the base of pa- 
hoehoe flows observed by (Keszthelyi, 1995), the formation 
of the "roller vortex" (Booth and Self, 1973), thermal erosion 
at low Reynolds numbers (Greeley et al., 1998), and the ob- 
servation, in lava channels, of magma temperatures greater 
than the eruption temperature at the vent (Kauahikaua et al., 
1998). Erosion by viscous heating in mantle convection was 
also investigated by Larsen et al. (1997) and Moore et al. 
(1998). 

Recently, the results of the studies of Polacci et al. (2001) 
suggest that the viscous heating may be responsible for the 
generation of a heterogeneous distribution of magma proper- 
ties inside the conduit. Moreover, the recent viscous gravity 
currents model of Vasilyev et al. (2001), based on the con- 
servation equations applied to a control volume, shows that 
the viscous dissipation exerts a strong influence in the vis- 
cous gravity currents. Finally, inadequacy of simple conduc- 
tive cooling models is shown by basal temperature measure- 
ments carried out by Keszthelyi (1995). Temperature mea- 
surements recorded at the base of many flows increases after 
some initial cooling. This fact can be easily understood on 
the light of the model presented here. 

We begin our study by investigating the limit case of flows 
in infinitely long tubes, and then we will describe the flow in 
tubes of finite length as being more representative of actual 
magma flows. We find that viscous heating effects cannot be 
generally neglected in typical magma flows and that the en- 
ergy and momentum equations cannot be decoupled. More- 
over, we find that the characteristic plug-like velocity profile 
observed in magma flows can be caused by viscous heating 
effects, even assuming Newtonian rheology. 

2 Model description 

In this paper, magma is assumed to be incompressible and 
approximated as an homogeneous fluid with constant den- 
sity, specific heat and thermal conductivity. The viscosity, 
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Fig. 1. Sketch of the inclined lava tube and prescribed boundary conditions. 



however, is temperature-dependent. We investigate the flow 
in a slab between two parallel boundaries separated by a dis- 
tance 2H and inclined of an angle a with respect to the hor- 
izontal (see Figure 1). For symmetry, the velocity is directed 
along the flow. At the tube walls (z = ±H) we impose a 
null velocity and a fixed temperature. This last assumption 
on the temperature is quite restrictive, and will be probably 
generalized in future works. However, for the case of tube 
of finite length, we show also results obtained by assuming 
thermal adiabatic condition at the walls, as a limit case. 

When the Peclet number is very high, the lubrication ap- 
proximation is valid: the transient term in the momentum 
equation is orders of magnitude smaller than the correspond- 
ing term in the energy equation (their ratio is of the order 
of the inverse of the Prandtl number). This means that the 
time scale for momentum relaxation is much shorter that the 
corresponding time scale for thermal relaxation, and the time 
evolution of the velocity and temperature profiles is controlled 
only by the energy equation. In these hypotheses, magma dy- 
namics in the tube are described by the following transport 
equations: 

-^+P9^na+—[n—\=0 (2) 

ox oz \ oz ) 

and 

c (^L + U—+W—\-k—+ (— V (3) 
^ \ dt dx dz ) dz 2 ^ \ dz ) 

where z is the coordinate transversal to the flow, x longitu- 
dinal coordinate, t time, (U, W) velocity of the fluid in the 
directions x and z respectively, p density, P pressure, p vis- 
cosity, g gravity acceleration, a slope, c specific heat at con- 
stant pressure, T temperature, and k thermal conductivity. 
The last term in (3) represents the heat generation by viscous 
dissipation. 

The characteristic length scales of this problem are the 
channel dimensions 2H (thickness) and L (length), the me- 
chanical relaxation length L m = UH 2 p/p , and the thermal 
relaxation length L t = UH 2 pc/k. For lava flows, typically 
L m /L <C 1, but L t /L 3> 1 and the approximation of in- 
finitely long tube is not valid. Only when L t /L -C 1, the 
approximation of infinitely long tube (by the thermal point 
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of view) is allowed. This limit case is simpler to handle, and 
it was widely studied in the past. For these reasons, we dis- 
cuss magma flow in infinitely long tubes, before studying the 
viscous heating effects in tubes of finite length. 

In silicate melts the dependence of the viscosity on the 
temperature is well described by the Arrhenius law (Shaw, 
1969; Danes, 1972; Park and Iversen, 1984; Baloga and Pieri, 
1986; Crisp and Baloga, 1994): 



/i = ha cxp(5/T) 



(4) 



where ha is a constant and B is the activation energy. In the 
interval of temperatures (T — To) /To <§C 1 (To is a reference 
temperature), eq. (4) is well approximated by the (Nahme's) 
exponential law: 



/i = p, exp[-/3(T - T )] 



(5) 



with (3 = B/T$ and p = p, A exp(B/T ). In the follow- 
ing, for simplicity, we adopt the exponential law, more suited 
for the analytical manipulation of the equations (Shaw, 1969; 
Dragoni, 1989; Costa and Macedonio, 2002). 

2. 1 Infinitely long pipe 

In this section we study the flow in an infinitely long chan- 
nel (slab) driven only by the component of the gravity force 
acting along the flow (all gradients are null along the flow). 
We assume no-slip conditions (U = 0), constant tempera- 
ture (T = T w ) at the walls (z = ±iT) and, for simplicity, 
we choose the reference temperature T = T w . Under these 
hypotheses, eq. (2) and (3) reduce to: 
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(6) 



(7) 



In a previous work, Pearson (1977) considered a Poiseuille 
flow between two horizontal parallel planes with isothermal 
boundaries (the temperature at the boundaries is the same as 
that of the injected fluid), negligible transversal gradients, 
and a given flow rate. Previously, a similar problem (the 
Hagen-Poiseuille flow) was numerically solved using ana- 
logue computers by Gruntfest et al. (1964), and later Eckert 
and Faghri (1986) studied the Couette flow. Like Pearson 
(1977), Gruntfest et al. (1964) found that during the tran- 
sient phase, in contrast to the steady-state, the temperature 
profile shows higher values near the walls than in the chan- 
nel centre. Moreover, Gruntfest et al. (1964) found that the 
solutions of the equations that they consider do not diverge 
only when the non-dimensional pressure gradient Q p (Q p = 
/3(dP/dx) 2 H 4 /kpo) is smaller than a critical value. In fact 
for Q p greater than this value, the temperature increases in- 
definitely in a finite time. Gruntfest et al. (1964) suggest that 
the behaviour of this model can help in understanding the ori- 
gin of turbulence (however this feedback could not produce 
an indefinite temperature increase in 2-D and 3-D systems). 



Next we investigate the time-dependent problem of a fluid 
with initial temperature T, and wall temperature T w = T , 
with T > T). We study the conditions for the existence 
of the solution, and the evolution to the steady-state solution 
(when it exists). After combining eq. (6) and (7), we rewrite 
the energy equation in the non-dimensional form: 



e = o 
e = B 

where 

= 0(T-T O ); C = 
f3(pg sinqfff 4 

6 = z ; 



Vr > 0, 
r = 0, 



C = ±i 
l < c < i 



t = kt/{pcH 2 ) 



B = m - To) 



(8) 



(9) 



(10) 



The parameter Q represents the non-dimensional "shear stress", 
and is known as the Nahme number based on the character- 
istic stress (instead of the velocity). 

From the results of Gruntfest et al. (1964) we know that 
when B = and Q > Q cr i t (Gcrit ~ 5.64 for slab flows), 
the solution of eq. (8) diverges in a finite time. In the case 
of Q < Gcrit, the solution of the transient phase, with two 
maxima near the walls, evolves for t w 1 to the steady-state 
solution with only one maximum in the channel centre. In 
Figure 3 and Figure 2 we report the temporal evolution of 
9 found numerically for the subcritical condition (Q = 5), 
and for the supercritical condition (Q = 8) respectively, both 
with B = 0. As shown in Figure 3, when Q < Q c hu f° r 
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Fig. 2. Temporal evolution of the non-dimensional temperature profile in 
the supercritical condition Q = 8 and B = 0. 

small t the temperature increases near the boundaries and 
for r w 1 its profile tends to the steady-state solution. In 
Figure 2, in the supercritical case, the temperature profile, 
shows two peaks near the walls that continuously increase, 
until they diverge when r reaches a critical value T cr u. 





Fig. 3. Temporal evolution of the non-dimensional temperature profile in 
the subcritical condition 5 = 5 and B = 0. 



The general case, with the initial fluid temperature greater 
than the temperature at the walls (B > 0), was previously 
studied by Fujita (1969) from a theoretical point of view, but 
is not yet analytically or numerically solved (to our knowl- 
edge). Fujita (1969) clarified the relationship between the 
transient and the steady-state phases with some theorems for 
the steady-state boundary value problem: 



vV 



e u = (x e 0) 



and for the corresponding time-dependent problem: 



— = V 2 U + e v 
at 



(t > o, x e ft) 



(ID 



(12) 



with initial condition u(t = 0) = a(x) (with a(x) continu- 
ous), and boundary condition u(d£Y) = 0. One of the theo- 
rems states that if eq. (1 1) has no solution, then the solutions 
of eq. (12) diverge in a finite time, or diverge for t — > oo. 
Another theorem states that if eq. (1 1) has more than one so- 
lution, 4> is the "smaller" solution, and ip is another solution 
not equal to </>, then: 

1. If a > tp, then the solution u of eq. (12) diverges, in a 
finite time, or for t — > oo. 

2. If a < ip, then the solution u of eq. (12) uniformly con- 
verges to <f), for t — > oo. 

From our numerical solutions of eq. (8), we observe that when 
B > 0, the range of values of Q below which the solution of 
eq. (8) converges is smaller than the corresponding range of 
Q for B = 0, but the profiles evolution is similar to the case 
with B = 0. We call Q cr it(B) the critical value of Q above 
which the solution diverges when the initial condition is B. 
An example is reported in Figure 5 for Q = 4 and B = 3 
where the temperature diverges in a finite time. Figure 6 
shows the relation between B and Q C rit(B) found numeri- 
cally. For B > we study the following two cases: in the 




Fig. 4. Temporal evolution of the non-dimensional temperature for Q = 2 
and B = 3. 
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Fig. 5. Temporal evolution of the non-dimensional temperature profile in 
the supercritical condition Q = 4 and B = 3. 



subcritical case (Q < Q cr it) we investigate the effect of the 
viscous heating on the temporal evolution of the temperature 
profile, and we compare the results to the case without heat 
generation (Q — 0); in the supercritical case, we investigate 
the relation between the critical time corresponding to the 
thermal blow-up T crit and Q. When the viscous heating is ne- 
glected, eq. (8) reduces to the classical heat equation. In fact 
from an initial condition with uniform temperature greater 
than that of the boundaries, the fluid starts cooling near the 
walls, until the wall temperature is reached in a time scale of 
the order of the characteristic time, i.e r ~ 1. By accounting 
for the viscous heating, as shown, for example, in Figure 4, 
the temporal evolution of the temperature can be described in 
four phases. The initial phase is characterized by the cooling 
near the walls; in the second phase the viscous heating near 
the walls produces the increase of the temperature above its 
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Fig. 6. Relation between B and Q cr it{B). 



initial value and the temperature profile assumes two max- 
ima in the external part of the profile and a minimum in the 
central part. During the third phase the two maxima migrate 
towards the central part of the channel, and the profile is char- 
acterized by a maximum temperature in the centre, that could 
be also greater than its initial value. In the fourth phase the 
fluid cools until the steady state solution is reached, in agree- 
ment with the theoretical results of Fujita (1969), previously 
reported. 

In the supercritical case (Q > Qcrit), we find an non- 
dimensional critical time T crit above which the solution di- 
verges. In this case, the temporal evolution of the tempera- 
ture profile is qualitatively similar to the corresponding case 
with B = 0, but with shorter T crit . The temperature increases 
near the walls until it diverges as r — > T CT it (similar to the 
case shown in Figure 5). As Q increases, T cr i t gets shorter. 
As an example, Figure 7 shows the relation between Q and 

Tcrit, for B = 0. 

From the physical point of view, when viscous heating is 
relevant, a transversal temperature gradient develops near the 
walls to dissipate the heat through the boundaries. This tem- 
perature increase produces a viscosity decrease with a con- 
sequent increase of the flow velocity. This velocity increase 
produces an increase in the velocity gradients near the walls, 
and a further local temperature increase. When Q > Q C r%t, 
the temperature increases without bound this feedback can- 
not converge and the steady-state solution does not exist. Of 
course, this has no physical meaning, and is related to the 
loss of validity of some assumptions such as the hypothe- 
sis of one-dimensional laminar flow. When the infinitely 
long pipe hypothesis is valid, according to Gruntfest et al. 
(1964), it is legitimate to assume that above the critical val- 
ues (G C rit, T cr it), other diffusive processes are triggered, al- 
lowing a more efficient heat and momentum transfer. 
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Fig. 7. Relation between r crit and Q for B = 0. 
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Instead, when Q < Q cr it, the temperature and velocity 
profiles tend to the steady-state solutions described by: 



dc, 2 



+ Qe\ 2 = 



(13) 



Eq. (13) was widely studied in the past, and it is well known 
that it has no solutions for Q > Q cr it, where Q cr it is a critical 
value (Joseph, 1964). Below the critical value, for each Q 
there are two solutions corresponding to different tempera- 
ture profiles: one at low temperature, and one at high tem- 
perature (see Figure 8). For some years, the existence of 
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Fig. 8. Non-Dimensional "shear-stress" Q vs. 9 max- (8 max is the maxi- 
mum temperature in the channel centre, at steady state.) 

multiple solutions corresponding to a given "shear-stress" Q 
was discussed in the literature, and the stability of the solu- 
tion in the higher branch was investigated. This depends on 
which is the controlling variable of the problem: the veloc- 
ity or the shear stress (John and Narayanan, 1997). More- 
over, the relevance of the boundary conditions to the multi- 
ple steady-states and the stability of the solution was investi- 
gated in the contest of asthenospheric shear flow (Yuen and 
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Schubert, 1977). A complete review on this argument un- 
til 1974 is reported by Sukanek and Laurence (1974), and 
more recently by Becker and McKinley (2000). In summary, 
in the Poiseuille flow between two parallel planes, when the 
non-dimensional "shear-stress" parameter Q is greater than 
a critical value Q crit (~ 5.64), eq. (13) has no solutions; for 
G = Qcrit it has one solution; and for Q < Gcrit it has two 
solutions, one of which (the one with greater temperature) 
may be unstable. For Q < G c hu typical profiles of the non- 
dimensional temperature 9, and velocity u/Uo that satisfy 
(6), (7) and the imposed boundary conditions, are character- 
ized by a maximum in the centre of the tube. 

2.2 Finite length tube 

From the thermal point of view, magma flows at high Peclet 
number cannot be described by assuming infinitely long tubes. 
At high Peclet numbers, the leading term in the left side 
of eq. (3) is the advective term, containing the longitudinal 
gradient of the temperature (e.g. Lawal and Kalyon, 1997). 
Including only the relevant terms, the transport equations 
(mass, momentum and energy balance) for finite length tubes, 
in non-dimensional form are: 



ud( = q (mass) 
du /dp . \ * a 

— = — — Sin a C e (momentum) 

n , do d 2 e AT , (du\ 2 _ e 

Pe u trrw +Na {dc) e (energy) 



(14) 



where?i = U/U* with U* = pgH 2 /p, , and p non-dimensional 
pressure (p = P/pgH); moreover: 



q = 



pgH 3 



Pe* = 



p 2 cgH 3 
kp, 



Pe 

q 



Na* = 



(3p 2 g 2 H 4 = Nao 
kjjL q 2 



(15) 



(16) 



(17) 



where Na* and Pe* are the Nahme and the Peclet num- 
bers based on the characteristic velocity U* respectively and 
Q = UH. The equation describing mass conservation is 
written in integral form. To solve the equations in (14), we 
assume no-slip conditions for the velocity and constant tem- 
perature T w at the walls. At the inlet we assume a parabolic 
velocity profile and constant temperature Tj n . Here, for sim- 
plicity and according to other authors, we consider the wall 
temperature as reference, i.e. T = T w and B = /3(T in —T w ). 
From an operative point of view, it would probably be more 
convenient to set the reference temperature equal to the tem- 
perature at the inlet. Of course, the physical behaviour of the 
flow is independent of the choice of the reference tempera- 
ture T ; this is obtained by rescaling all the parameters. As 



an example, by assuming T = T in , the viscosity needs to be 
rescaled by a factor e B , etc. 

Recently, Lawal and Kalyon (1997) found analytical so- 
lutions of the simpler problem with constant viscosity and 
showed typical temperature profiles with two peaks near the 
walls. In our case, since we consider temperature-dependent 
viscosity, we expect a feedback on the velocity profile. In 
fact, using asymptotic solution method for high Nahme num- 
bers, Ockendon (1979) solved the eq. (14) taking in to con- 
sideration the full advective term. These solutions describe 
the evolution of a flow from the Poiseuille regime with uni- 
form temperature at the conduit inlet to a plug regime char- 
acterized by two peaks in the temperature profile near the 
boundaries far from the inlet. Since the ratio L t /L m (equal 
to the Prandtl number) is very high, it is legitimate to assume 
that the flow is fully developed, by the mechanical point of 
view, almost at the tube inlet. The characteristic distance 
from the inlet where the viscous heating becomes relevant, 
is £* - PeNa- 3 ? 2 (Ockendon, 1979) (using data from Ta- 
ble 1 this gives 50-=-5000m for a conduit width of 5m). Under 
more general conditions, Schneider (1976) found numerical 
solutions of the viscous heating problem, showing clearly 
that the deviation from the Newtonian behaviour may have 
a thermal origin in fluids with temperature-dependent vis- 
cosity. When the viscous dissipation becomes important, its 
effect overcomes the thermal cooling from the walls produc- 
ing a maximum in the temperature profile near the walls and 
a consequent plug-like velocity distribution. These effects 
increase with distance from the inlet and, unexpectedly, the 
maximum temperature gets more pronounced when the tem- 
peratures at the walls get lower. Similar results were obtained 
by Galili et al. (1975) using perturbative methods: they found 
that the flow rate is not proportional to the longitudinal pres- 
sure gradient and the apparent viscosity of a Newtonian fluid 
depends on the shear rate because of the viscous dissipation. 

We solve eq. (14) by using a finite-difference method with 
an implicit scheme for the integration along direction £; the 
pressure gradient is iteratively adjusted at each step in order 
to satisfy mass conservation. For symmetry we investigate 
only half a channel (0 < C < !)• We find that the process is 
controlled by four parameters: a Peclet number Pe, a Nahme 
number Na*, the non-dimensional flow rate q, and the non- 
dimensional input temperature B. However, since we focus 
on magma flows in tubes, we set Pe = 10 7 as a typical value, 
and we perform a parametric study on the others. 

In Figure 9 and 10, we show a solution of eq. (14) for dif- 
ferent non-dimensional distance from the inlet: from £ = 
to £ = 10 4 . Starting with uniform temperature (6 = 0) 
and parabolic velocity profile at the inlet, the flow evolves 
gradually to a plug-like velocity profile with two symmetric 
peaks in the temperature distribution. In this case we used 
Na* = 100, q — 1 and B = and the viscous dissipation 
effects become important only for high values of £. Instead 
in Figure 11 and 12 the solutions for Na* = 1000, q = 1 
and B = are shown; these results are qualitatively similar 
to the previous ones but the temperature peaks are more pro- 
nounced and the length scale for the development of the plug 
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flow is lower than before (practically, the transition occurs 
for£ = 1000). 

When the viscous heating is important with B > 0, the 
typical temperature profile shows low values at the walls with 
peaks near the walls, as reported in Figure 13. As shown in 
Figure 14, starting with a parabolic velocity profile at the 
inlet (£ = 0), the flow migrates to a plug-like regime down- 
stream. By increasing B, the peak in the temperature profile 
moves towards the centre of the channel. 

In general, viscous heating becomes important when ei- 
ther the Nahme number (Na*) or the non-dimensional flow 
rate (q) increase. For comparison, in Figure 15 and 16, we 
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Fig. 9. Longitudinal evolution of the non-dimensional temperature profile 
Pe = 10 7 , Na* = 100, q = 1, B = for different non-dimensional 
distance from the inlet (10~ 3 § = 0, 2, 4, 6, 8, 10). 




Fig. 10. Longitudinal evolution of the non-dimensional velocity profile 
Pe = 10 7 , Na* = 100, q = 1, B = for different non-dimensional 
distance from the inlet (10~ 3 § = 0, 2, 4, 6, 8, 10). 

show the evolution of the temperature and velocity profiles in 
a case similar to Figure 9 and 10, but with adiabatic thermal 
conditions at the walls. In this case, as expected, the temper- 
ature increases more than the case of isothermal walls, and 
the the velocity profile shows a more pronounced plug-like 
behaviour. 




Fig. 11. Longitudinal evolution of the non-dimensional temperature profile 
Pe = 10 7 , Na* = 1000, q = 1, B = for different non-dimensional 
distance from the inlet (£ = 0, 200, 400, 600, 800, 1000). 
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Fig. 12. Longitudinal evolution of the non-dimensional velocity profile 
Pe = 10 7 , Na* = 1000, q = 1, B = for different non-dimensional 
distance from the inlet (£ = 0, 200, 400, 600, 800, 1000). 
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Fig. 13. Longitudinal evolution of the non-dimensional temperature profile 
Pe = 10 7 , Na* = 1000, q = 4, B = 4 for different non-dimensional 
distance from the inlet (10~ 3 £ = 0, 2, 4, 6, 8, 10). 
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Fig. 14. Longitudinal evolution of the non-dimensional velocity profile 
Pe = 10 7 , Na* = 1000, c? = 4, B = 4 for different non-dimensional 
distance from the inlet (10~ 3 £ = 0, 2, 4, 6, 8, 10). 
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Fig. 15. Longitudinal evolution of the non-dimensional temperature pro- 
file Pe = 10 7 , Na* = 100, q = 1, adiabatic walls, for different non- 
dimensional distance from the inlet (10~ 3 £ = 0, 2, 4, 6, 8, 10). 



2.3 Velocity field stability and secondary flows 

The energy equation of fluids with temperature-dependent 
viscosity in infinitely long flows is similar to the equation 
which governs chemical explosive processes in some mate- 
rials (described by the Frank-Kamenetski equation, [1939]). 
In this last case, when the parameter corresponding to Q is 
greater than its critical value, ignition occurs. In fluids, how- 
ever, other multidimensional degrees of freedom can be ac- 
tivated without the occurrence of extreme events. For ex- 
ample, in flows between two parallel planes, cylindrical sec- 
ondary flows can develop near the walls, or toroidal sec- 
ondary flows can occur in circular pipes. 

Due to the strong coupling between viscosity and temper- 
ature, the thermal instability generated by viscous heating 
may trigger an instability in the velocity field, which can- 
not be predicted by a simple isothermal newtonian models. 
When the viscous heating produces sharp peak in the tem- 
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Fig. 16. Longitudinal evolution of the non-dimensional velocity profile 
Pe = 10 7 , Na* = 100, q = 1, adiabatic walls, for different non- 
dimensional distance from the inlet (10~ 3 £ = 0, 2, 4, 6, 8, 10). 



perature profile near the walls, with consequent strong in- 
crease in the viscosity gradient, a triggering of instabilities 
and the transition to secondary flows (with more efficient 
thermal and mechanical diffusion) is possible. The stability 
of the plane Couette flow was recently re-examined by Yueh 
and Weng (1996), who improve the results previously found 
by Sukanek et al. (1973). The flow shows two different in- 
stability modes: the former arises in the non-viscous limit, 
and the latter is due to the viscosity stratification (Yueh and 
Weng, 1996). For this last instability mode, it was numeri- 
cally demonstrated that the critical Reynolds number (Re c ), 
above which the flow becomes turbulent, decreases with the 
increase of the Nahme number (Na), that is with the viscous 
heating. This behaviour is confirmed by recent experiments 
performed by White and Muller (2000). In these experi- 
ments, the authors use a temperature-dependent fluid (glyc- 
erin) and a Taylor-Couette device which allows the tracking 
of vortex by a laser particle tracer. Results clearly show that, 
above a critical Nahme number, an instability appears at a 
Reynolds number one order of magnitude lower than the cor- 
responding Reynolds number predicted for isothermal flow. 

Moreover, the triggering of the instability in the laminar 
flow and the activation of new secondary flows are confirmed 
by our 2-D direct numerical simulations of the complete Navier- 
Stokes equations, which will be the subject of a further work. 



3 Implications for lava flows 

In this section we apply the previous theoretical results to the 
study of lava and magma flows in tubes. We take the typical 
magma parameters reported in Table 1 as representative, but 
the obtained results can be easily generalized. 

Following Bruce and Huppert (1989), T w is the tempera- 
ture at the wall, which is defined as the temperature at which 
crystallization of magma ceases the flow. Here, we estimate 
T w as the temperature in the mid-range between T in and the 
solidification temperature T s : therefore, we use Tq = T w 
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Table 1. Parameters characteristic of lava flows. 
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1273K. From this value of T w and from the values of (3 and 
Ti n , reported in Table 1, we have B — (3(Ti n — T w ) w 4 
and p — ^{Tin) • e B « 30, 000 Pas. Using the other values 
reported in Table 1, we obtain the following characteristic 
non-dimensional numbers: 



q = fiaQ/pgH 3 
Na = p U 2 (3/k 
Pe = pcUH/k 
Pr Q = cpo/k 
Na - p(T m )U 2 f3/k 
Re = P UH/n{T in ) 
Pr = cn{T in )/k 



1 -j- 10 

6 • 10 2 6 • 10 4 

3-10 5 -^3-10 7 

1.5 • 10 7 

1 + 10 3 

1 -j- 100 

3 • 10 5 



(18) 



The last three numbers in eq. (18) are the Nahme, the Reynolds 
and the Prandtl numbers based on /j,(T in ), respectively, and 
are reported just for completeness. 

The values of the non-dimensional numbers reported in 
(18) allow for the application of the model described above 
for lava and magma flows. In fact, the high Peclet number 
results in a very fast advective transport and validates the as- 
sumption of the lubrication theory. Concerning the effects 
of viscous heating on magma flows, we found that for the 
lower value of the non-dimensional flow rate considered in 
our study (q = 1) and Na* = Na /q 2 < 1000, the cooling 
effects prevail on the viscous heating. This becomes domi- 
nant for Na* <~ 10 4 showing velocity and temperature pro- 
files similar to those reported in Figure 13 and 14. Instead, 
at high flow rates (q = 10), the corresponding values of Na* 
are lower: the viscous heating effects are dominant even at 
Na* r~j 100 whereas for Na* < 10 cooling effects prevail. 
Since these conditions are usual in magma flows, it should 
be possible to observe viscous heating effects in the natural 
environment. 

For magma flows in eruptive conduits, recently, Polacci 
et al. (2001) emphasized the importance of the viscous heat- 
ing in magma flow during volcanic eruptions to explain the 
creation and the discharge of two different varieties of pumice. 
In their scheme, one type of pumice originates in the region 
with greater temperature and higher strain rate near the con- 
duit walls, whereas the other type is generated in the central 
part of the conduit with lower strain rate and temperature. 

For lava flows, field evidences corresponding to the start- 
ing of secondary flows, previously described, are possibly 
represented by the "roller vortex" phenomenon (Booth and 



Self, 1973) and by the thermal erosion observed by Gree- 
ley et al. (1998) where the authors find "unequivocal evi- 
dence for thermal erosion" in lava tubes at Cave Basalt, Mt. 
St.Helens for which the dynamic analysis indicates laminar 
flows. It is known that the turbulent flow is more effective 
in erosion than the laminar flow because of the higher heat 
transfer rate (Hulme, 1973; Greeley et al., 1998); moreover 
preliminary studies show that purely thermal erosion by lam- 
inar flow is very difficult unless the substrate is of much 
lower melting temperature than the eroding fluid (Greeley 
et al., 1998). Field studies report that during the Kilauea 
eruption in 1994, lava eroded 5 meters of the basaltic sub- 
strate with an average erosion rate of 10 cm/day at low Reynolds 
numbers between 16 and 64 (Kauahikaua et al., 1998). Such 
high erosion rates highlight the difficulties that arise when 
trying to explain thermal erosion in laminar flows. More- 
over, Kauahikaua et al. (1998) show that the Jeffreys equa- 
tion that relates the velocity and the flow depth in a lami- 
nar flow is less adequate than other relations valid for turbu- 
lent flows such as the equation of Goncharov-Chezy, even for 
flows at low Reynolds numbers. Concerning the temperature, 
Kauahikaua et al. (1998) report the interesting and puzzling 
observations of some temperature measurements greater than 
the magma temperature at the vent and even greater that the 
upper limit of the temperature range of the used radiometer 
(about 1200 °C). An explanation of this phenomenon could 
be the local increase of temperature near the walls caused by 
the viscous heating, as previously described. Other authors 
also report the evidence of vortices associated with the begin- 
ning of turbulence in active lava flows (Keszthely and Self, 
1998) and the onset of complex flow patterns (Greeley, 1987; 
Lipman and Banks, 1987; Crisp and Baloga, 1994; Calvari 
and Pinkerton, 1998). The inadequacy of pure laminar flows 
to explain field observation such as mixing and streamlines 
breaking has been well known since the sixties when the not 
widely accepted term of "disrupted flow" was introduced to 
describe a flow with characteristics between pure laminar and 
turbulent (for a discussion see Baloga et al., 1995). Moreover 
basal temperature measurements carried out by Keszthelyi 
(1995) show that temperature at the base of some flows in- 
creases after some initial cooling. This behaviour cannot be 
predicted by a simple conductive model but can be explained 
using a correct physical model which correctly describes vis- 
cous heating effects. 

At Etna, during the eruption of 1971, secondary rotational 
flows were observed near the walls in flows confined between 
levees or in deep channels (Booth and Self, 1973). These sec- 
ondary flows consist of two elicoidal patterns symmetrical to 
the centre of the flow. Booth and Self (1973), because of 
the low Reynolds number, exclude any turbulent flow, and 
try to explain these rotational flows invoking other causes 
and other, sometimes weak, argumentations. We think that 
the observed phenomena reported by Booth and Self (1973) 
could be interpreted as the secondary flows triggered by a 
pronounced viscous heating. The rotational flows, which al- 
low the mixing of the fluid along sections orthogonal to the 
direction of the flow, invoked in the model of lava flows with 
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two thermal components by Crisp and Baloga (1990), could 
have a similar interpretation. 

Moreover, exploration and investigations on terrestrial plan- 
ets and satellites demonstrated that their lava channels are 
typically larger than the corresponding channels on the Earth. 
To account for the large sizes of lunar and martian chan- 
nels Hulme (1973, 1982) and Carr (1974) proposed that lavas 
eroded the ground over which they flowed (see also Baloga 
et al. (1995)). Our model could be used to investigate whether 
the action of secondary flows caused by viscous heating could 
have eroded the preflow surface more efficiently than laminar 
flows. 

We wish to remark that the proposed model is applicable 
to slab flows, but it easily generalizable to flows in circular 
pipes. Moreover, modeling of magma and lava flows needs 
also the study of open channels with free surface and differ- 
ent types of thermal boundary conditions, not considered in 
the present study. In particular, a more realistic description 
should account for the temperature variations of the ambi- 
ent medium, that needs the introduction of additional con- 
trol parameters such as the Nusselt number which measures 
how much the heat transferred through the boundaries is con- 
ducted away. 

Finally, an accurate study of the flow instabilities and the 
triggering of secondary flows needs the solution of the com- 
plete transport equations in 2-D or 3-D, which is the subject 
of a further work in preparation. 

4 Conclusions 

The effects of viscous heat generation in fluids with strongly 
temperature-dependent viscosity such as silicate melts are in- 
vestigated. These effects can play an important role in the 
dynamics of magma flow in conduits and lava flows in tubes. 
In fact, viscous friction generates a local increase in temper- 
ature near the tube walls, with consequent increase of the 
fluid velocity because of the viscosity decrease. In the case 
of infinitely long tubes, one-dimensional models predict that 
viscous heating can lead to a positive feed-back known as 
"thermal runaway". Actually, as described by 2-D and 3- 
D models, this feed-back causes the activation of other de- 
grees of freedom, with the production of local instabilities 
and the triggering of secondary flows. However, the typical 
high value of the Peclet number in magma flows does not al- 
low the assumption of an infinitely long tube (from the ther- 
mal point of view), and a model for magma flow in tubes of 
finite length is needed, as described in the paper. In general, 
this process is controlled by the value of the Nahme number, 
the flow rate, the Peclet number, and the thermal conditions 
at the inlet. By adopting such a model, and typical values of 
the parameters for lava flows in tubes, we find that starting 
from a constant temperature and parabolic velocity profile at 
the tube inlet, viscous heating causes the increase of the tem- 
perature near the walls, with a consequent local viscosity de- 
crease. This can lead to the formation of a plug-like velocity 
profile, commonly observed in lava flows, and explained only 
by assuming a Binghamian rheology. Moreover, the presence 
of an inflex in the velocity profile near the walls can lead to 
the formation of local flow instabilities, even at low Reynolds 
numbers. The assumption of isothermal boundary conditions 
at the walls, although simplifying, is probably too restrictive; 



however, by assuming adiabatic conditions, the qualitative 
character of the flow does not seem to change dramatically 
in tubes of finite length. We hope that this work will provide 
a motivation for further theoretical and field investigations 
of viscous heating effects in magma flows. In particular, the 
formation of vortices in active lava flows at low Reynolds 
numbers, thermal erosion, temperature and velocity profile. 
Moreover, the effects viscous dissipation in volcanic con- 
duits could have an important role on the dynamics of the 
both effusive and explosive eruptions. 
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